﻿using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;
using System.Threading.Tasks;

namespace AnalysisSystem
{
    class ArpsMethod
    {
        //注意所有的循环变量开始值及结束值
        public static bool GoldArps(int doti, int dote, List<double> qo, List<double> tx, double ArpsN, double intercept, double slope)
        {
            //doti = 起始点 dote = 结束点
            //ln(1/Q^n) = ln(1/Q^n) + (Di*n/Di^n)*(t-ti)
            double intervalA, intervalB, intercept1, intercept2;
            double y1 = 0, y2 = 0, y3, slope1, slope2, N1, N2, R = 0;
            int dot = dote - doti + 1;
            double[] QN1 = new double[dot];
            double[] QN2 = new double[dot];
            double[] tn = new double[dot];

            //寻找最佳n
            intervalA = 0;
            intervalB = 1;
            while(Math.Abs(intervalB - intervalA) >= 0.00001)
            {
                N2 = intervalA + 0.618034 * (intervalB - intervalA);
                N1 = intervalA + 0.381966 * (intervalB - intervalA);
                int j = 0;
                for(int i=doti; i<dote; i++, j++)
                {
                    tn[j] = tx[i] - tx[doti];
                    QN2[j] = 1 / (Math.Pow(qo[i], N2));
                    QN1[j] = 1 / (Math.Pow(qo[i], N1));
                }
                Regress1(j, tn, QN1, out slope1, out intercept1);
                Regress1(j, tn, QN2, out slope2, out intercept2);
                for(int i=0; i<j; i++)
                {
                    y1 = y1 + Math.Pow((QN1[i] - intercept1 - slope1 * tn[i]), 2);
                    y2 = y2 + Math.Pow((QN2[i] - intercept2 - slope2 * tn[i]), 2);
                }
                if(y2 > y1)
                {
                    intervalB = N2;
                    N2 = N1;
                    N1 = intervalA + 0.381966 * (intervalB - intervalA);
                }    
                else
                {
                    intervalA = N1;
                    N1 = N2;
                    N2 = intervalA + 0.618034 * (intervalB - intervalA);
                }
            }
            ArpsN = 0.5 * (intervalB + intervalA);
            y3 = 0.5 * (y1 + y2);

            //与特殊的指数递减作比较
            int j = 0;
            for(int i=doti; i<dote; i++)
            {
                if(qo[i] > 0)
                {
                    tn[j] = tx[i] - tx[doti];
                    QN1[j] = Math.Log(qo[i]);
                    j++;
                }
            }
            Regress1(j, tn, QN1, out slope1, out intercept1);
            y1 = 0;
            for(int i=0; i<j; i++)
            {
                y1 = y1 + Math.Pow((QN1[i] - intercept1 - slope1 * tn[i]), 2);
            }
            if(y3 > y1)
            {
                y3 = y1;
                ArpsN = 0;
            }

            //与特殊的调和递减比较
            j = 0;
            for(int i=doti; i<dote; i++)
            {
                if(qo[i] > 0)
                {
                    tn[j] = tx[i] - tx[doti];
                    QN2[j] = 1 / qo[i];
                    j++;
                }
            }
            Regress1(j, tn, QN2, out slope2, out intercept2);
            y2 = 0;
            for(int i=0; i<j; i++)
            {
                y2 = y2 + Math.Pow(QN2[i] - intercept2 - slope2 * tn[i], 2);
            }
            if(y3 > y2)
            {
                y3 = y2;
                ArpsN = 1;
            }

            //根据最佳n重新计算理论Q(t)
            if(ArpsN < 0.00001)
            {
                j = 0;
                for(int i=doti; i<dote; i++)
                {
                    if(qo[i] > 0)
                    {
                        QN1[j] = Math.Log(qo[i]);
                        j++;
                    }
                }
            }
            else
            {
                j = 0;
                for(int i=doti; i<dote; i++)
                {
                    if(qo[i] > 0)
                    {
                        QN1[j] = 1 / Math.Pow(qo[i], ArpsN);
                    }
                }
            }
            Regress1(j, tn, QN1, out slope, out intercept, R);
            return true;
        }

        private static bool Regress1(int n, double[] x, double[] y, out double b, out double a, double R = 0)
        {
            //y(x) = A+B*x
            //R=Relative Coeff
            if (n == 0)
            {
                b = 0; a = 0;
                return false;
            }
            double sx = 0, sxx = 0, sy = 0, syy = 0, sxy = 0, u = 0, v = 0;
            for(int i=0; i<n; i++)
            {
                sx += x[i];
                sy += y[i];
                sxx += x[i] * x[i];
                syy += y[i] * y[i];
                sxy += x[i] * y[i];
            }
            u = (sxx - sx * sx / n) + 0.000001;
            b = (sxy - sx * sy / n) / u;
            a = (sy - b * sx) / n;
            v = Math.Abs(syy - sy * sy / n) + 0.000001;
            R = Math.Abs(sxy - sx * sy / n) / Math.Sqrt(Math.Abs(u * v));
            return true;
        }

    }
}
